#Uses devel version of AlphaSimR '0.3.0'
# > devtools::install_bitbucket("hickeyjohnteam/AlphaSimR", args = '--library=.')
library(AlphaSimR)
#Create MaCS command for 10 breed split
#'size' should be 2 time the number of individuals and a factor of 10

# options(scipen=999)
args = commandArgs(trailingOnly=TRUE)
# args = c("10", "100", "1")

nChr = as.numeric(args[1])
nSnp = as.numeric(args[2])
maplength = as.numeric(args[3])

print(paste(nChr, nSnp, maplength))

nInd = 100
nSnpPerChr =  rep(nSnp, nChr) #Must be larger than nQtl

# founderPop = runMacs(nInd=nInd, nChr=nChr, segSites=nSnpPerChr, species = "CATTLE")
founderPop = runMacs(nInd=nInd, nChr=nChr, segSites=nSnpPerChr, inbred=TRUE, species="WHEAT")

# ---- Set Simulation Parameters ----

SP = SimParam$new(founderPop)
SP$addTraitA(nQtlPerChr=nSnpPerChr, mean=0, var=1)
SP$addSnpChip(nSnpPerChr=nSnpPerChr)
SP$setVarE(h2=.4)

# SP$setRecRatio(maplength)
tmp = lapply(SP$genMap,function(x) (maplength*x)) # Modifies the recombination rate based on the paternal chromosome.
SP$switchGenMap(tmp)


# genRef = newPop(founderPop)
# genFinal = randCross(genRef, nCrosses = nInd, nProgeny = 1) 

genRef = newPop(founderPop)
gen1half = selectCross(genRef, nInd = nInd, nCrosses = nInd, nProgeny = 1)
genFinal = makeDH(gen1half, nDH = 1)


write.table2 = function(x, file, ...) {
    write.table(x, file, col.names=F, row.names=F,quote=F,...)
}

writeGenotypes = function(subPop, suffix){
    genotypes = AlphaSimR::pullSnpGeno(subPop, 1)
    haplotypes = AlphaSimR::pullSnpHaplo(subPop, 1)
    pedigree = data.frame(subPop@id, subPop@father, subPop@mother, stringsAsFactors = FALSE)

    write.table2(cbind(as.numeric(pedigree[,1]), genotypes), paste0("genotypes.", suffix))
    write.table2(cbind(as.numeric(rep(pedigree[,1], each = 2)), haplotypes), paste0("haplotypes.", suffix))
}
writeGenotypes(genFinal, "children")
writeGenotypes(genRef, "parents")

output = c(genRef, genFinal)
pedigree = data.frame(output@id, output@father, output@mother, stringsAsFactors = FALSE)
write.table2(pedigree, "pedigree.txt")

map = data.frame(rep(1:nChr, each = nSnp), 0, rep(1:nSnp, times = nChr))
write.table2(map, "map.txt")

training_half = selectCross(genRef, nInd = 100, nCrosses = 1000, nProgeny = 1)
training = makeDH(training_half, nDH = 1)

write.table2(cbind(training@id, AlphaSimR::pullSnpGeno(training, 1)), "training_genotypes.txt")

